home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / simplx.arc / SIMPLEX.PAS < prev   
Pascal/Delphi Source File  |  1985-09-27  |  10KB  |  381 lines
  1. {
  2.     This routine performs least squares curve fitting to arbitrary
  3.    functions using the simplex method. This algorithm was derived
  4.    from an article in the May 1984 issue of BYTE magazine written
  5.    by Marco S. Caceci and William P. Cacheris, Florida State Univ.
  6.      Refer to that article for an indepth discussion on the mech-
  7.    anics of the program. Also see Nelder J.A. and R. Mead,
  8.    Computer Jou. 7, 308 (1965) and L.A. Yarbro and S.N. Deming,
  9.    Anal. Chim. Acta 74, 391 (1974).
  10.  
  11.                                    Writen for Turbo Pascal JUNE 1985
  12.                                    by James G. Landowski 72167,2736
  13.                                       Canoga Park, CA.
  14. }
  15.  
  16. {
  17.      Changed 2 to nvpp in several places.
  18.      Corrected relative error to allow for negative values
  19.  
  20.      Lew Paper
  21.      9/27/85
  22. }
  23.  
  24. const m     = 2;        { no. of parameters in equation to fit         }
  25.       nvpp  = 2;        { no. of variables per data point to solve for }
  26.       n     = 3;        { points in simplex ( m + 1 )                  }
  27.       mnp   = 200;      { max. no. of data points that can be fit      }
  28.       alfa  = 1.0;      { reflection coeff.   > 0    }
  29.       beta  = 0.5;      { contraction coeff.  0 => 1 }
  30.       gamma = 2.0;      { expansion coeff.    > 1    }
  31.       root2 = 1.414214;
  32.  
  33. type  vector  = array[1..n] of real;
  34.       datavec = array[1..nvpp] of real;
  35.       index   = 0..255;
  36.  
  37. var   done    :boolean;
  38.       i,j     :index;
  39.       hi,lo   :array[1..n] of index;
  40.       np,
  41.       maxiter,
  42.       numiter : integer;
  43.       next,
  44.       center,
  45.       mean,
  46.       error,
  47.       maxerr,
  48.       p,q,
  49.       step    :vector;
  50.       simp    :array[1..n] of vector;
  51.       data    :array[1..mnp] of datavec;
  52.       fname   :string[14];
  53.       din,
  54.       dout    :text;
  55.       abs_max :REAL; {L.P.}
  56. { ----------------------------------------------------------------------}
  57. {  Define the function to be fitted, 'C' is the coeff. vector,
  58.      'X' is the independent variable vector. This equation is the
  59.      example used in the BYTE article.
  60. }
  61.  
  62. function fnc(C:vector; X:datavec) :real;
  63.   begin
  64.       fnc := C[1] * X[1] /( C[2] + X[1] )
  65.   end;
  66.  
  67. { --------------------------------------------------------------------- }
  68. procedure sumresiduals (var s :vector);
  69.  
  70. { computes the sum of squares of the residuals }
  71.  
  72. var i:index;
  73.  
  74. begin
  75.   s[n] := 0.0;
  76.   for i := 1 to np do       { do for all data points }
  77.   begin
  78.     s[n] := s[n] + sqr( fnc(s,data[i])-data[i,nvpp] )
  79.   end
  80. end;
  81.  
  82. { ----------------------------------------------------------------------- }
  83. procedure inputdata;
  84.  
  85. { gets all input data from a disk file. The filename is prompted for
  86.   at runtime. An example of input data (from BYTE article) is:
  87.  
  88.   100               maximum no. of iteration in search
  89.   .2 3              starting coordinates of simplex
  90.   .1 1              starting increments
  91.   1e-4 1e-4 1e-4    maximum error
  92.   1.68 .172         data ( 'mnp' is max. no. of entries )
  93.   3.33 .25           |
  94.   5.00 .286          |
  95.   6.67 .303          |
  96.   10.0 .334          |
  97.   20.0 .384          V
  98. }
  99.  
  100. begin
  101.    writeln(dout);
  102.    writeln(dout,' Data file used for this run: ',fname);
  103.    readln(din,maxiter);
  104.    writeln(dout,' Maximum no. of iterations: ',maxiter:5);
  105.  
  106.    writeln(dout);
  107.    writeln(dout,' Starting coordinates: ');
  108.    write(dout,'    ');
  109.    for i := 1 to m do
  110.      begin
  111.       read(din,simp[1,i]);
  112.       write(dout,simp[1,i],' ');
  113.      end;
  114.    writeln(dout);
  115.    readln(din);
  116.  
  117.    writeln(dout);
  118.    writeln(dout,' Starting steps:');
  119.    write(dout,'    ');
  120.    for i := 1 to m do
  121.      begin
  122.       read(din,step[i]);
  123.       write(dout,step[i],' ');
  124.      end;
  125.    writeln(dout);
  126.    readln(din);
  127.  
  128.    writeln(dout);
  129.    writeln(dout,' Maximum errors:');
  130.    write(dout,'    ');
  131.    for i := 1 to n do
  132.      begin
  133.       read(din,maxerr[i]);
  134.       write(dout,maxerr[i],' ');
  135.      end;
  136.    writeln(dout);
  137.    readln(din);
  138.  
  139.    writeln(dout);
  140.    writeln(dout,' Data:');
  141.    np := 0;
  142.    while NOT EOF(din) do
  143.      begin
  144.        np :=succ(np);
  145.        write(dout,'     #',np:2);
  146.          for j := 1 to nvpp do
  147.            begin
  148.              read(din,data[np,j]);
  149.              write(dout,data[np,j])
  150.            end;
  151.         readln(din);
  152.         writeln(dout);
  153.      end
  154.   end;
  155. { ----------------------------------------------------------------------- }
  156.  
  157. procedure outputdata;
  158.  
  159. var y, dy, sigma :real;
  160. {    pltout       :text;}
  161.  
  162. begin
  163.  
  164. {  Open output files for plot data when needed.   }
  165.  
  166. {   assign(pltout,'SIMPLEX.PLT'); }
  167. {   rewrite(pltout);              }
  168.  
  169.    writeln(dout);
  170.    writeln(dout,' Completion after ',numiter:5,' iterations');
  171.    writeln(dout);
  172.    writeln(dout,' Final Simplex:');
  173.    for j := 1 to n do
  174.      begin
  175.        write(dout,'    ');
  176.        for i := 1 to n do
  177.           write(dout,simp[j,i]:12,'  ');
  178.        writeln(dout);
  179.      end;
  180.    writeln(dout);
  181.    writeln(dout,' Mean(s): ');
  182.    for i := 1 to n do
  183.      begin
  184.       if i <> n then writeln(dout,' C',i:1,'        ',mean[i]);
  185.       if i = n  then writeln(dout,' Residual: ',mean[i]);
  186.      end;
  187.    writeln(dout);
  188.    writeln(dout,' Fractional Error:');
  189.    write(dout,'    ');
  190.    for i := 1 to n do
  191.       write(dout,error[i]);
  192.    writeln(dout);
  193.    writeln(dout);
  194.    sigma := 0.0;
  195.    for i := 1 to np do
  196.      begin
  197.        y := fnc(mean,data[i]);
  198.        dy := data[i,nvpp]-y;
  199.        sigma := sigma+sqr(dy);
  200.      end;
  201.    sigma := sqrt(sigma/np);
  202.    writeln(dout,' Std. Dev. of error: ',sigma);
  203.    sigma := sigma/sqrt(np-m);
  204.    writeln(dout,' Estimate of error : ',sigma);
  205.    writeln(dout);
  206.    writeln(dout,'  #       x          y            y"         dy');
  207.    for i := 1 to np do
  208.      begin
  209.        y := fnc(mean,data[i]);
  210.        dy := data[i,nvpp]-y;
  211.        writeln(dout  ,i:3,' ',data[i,1]:12,' ',data[i,nvpp]:12,' ',y:12,' ',
  212.                dy:12);
  213. {      writeln(pltout,i:3,' ',data[i,1]:12,' ',data[i,nvpp]:12,' ',y:12,' ',
  214.                dy:12); }
  215.      end;
  216.    writeln(dout);
  217.    writeln(dout,'  -Done- ');
  218. end;
  219.  
  220. { ----------------------------------------------------------------------- }
  221.  
  222. procedure first;
  223.   begin
  224.    writeln(dout);
  225.    writeln(dout,' Starting SIMPLEX');
  226.    writeln(dout);
  227.    for j := 1 to n do
  228.      begin
  229.       write(dout,' simp[',j:1,']  ');
  230.       for i := 1 to n do
  231.          write(dout,simp[j,i]:12,' ');
  232.       writeln(dout)
  233.      end;
  234.    writeln(dout);
  235.   end;
  236.  
  237. { ----------------------------------------------------------------------- }
  238.  
  239. procedure newvertex;
  240.  
  241. begin
  242. {   write(dout,' -> ',numiter:4); }
  243.    for i := 1 to n do
  244.      begin
  245.        simp[hi[n],i] := next[i];
  246. {      write(dout,next[i]) }
  247.      end;
  248. {   writeln(dout)          }
  249. end;
  250.  
  251. { ----------------------------------------------------------------------- }
  252.  
  253. procedure order;
  254.  
  255. var i,j :integer;
  256.  
  257. begin
  258.    for j := 1 to n do
  259.      begin
  260.       for i := 1 to n do
  261.         begin
  262.          if simp[i,j] < simp[lo[j],j] then lo[j] :=i;
  263.          if simp[i,j] > simp[hi[j],j] then hi[j] :=i;
  264.         end
  265.      end
  266. end;
  267.  
  268. { ----------------------------------------------------------------------- }
  269. { Main }
  270. begin
  271.  
  272.    write('Filename: ');
  273.    readln(fname);
  274.    assign(din,fname);
  275.    reset(din);
  276.    assign(dout,'simplex.out');
  277.    rewrite(dout);
  278.  
  279.    inputdata;
  280.  
  281.    sumresiduals(simp[1]);
  282.  
  283.    for i := 1 to m do
  284.      begin
  285.       p[i] := step[i]*(sqrt(n)+m-1)/(m*root2);
  286.       q[i] := step[i]*(sqrt(n)-1)/(m*root2)
  287.      end;
  288.    for i := 2 to n do
  289.      begin
  290.       for j := 1 to m do simp[i,j] := simp[1,j] + q[j];
  291.       simp[i,i-1] := simp[1,i-1]+p[i-1];
  292.       sumresiduals(simp[i])
  293.      end;
  294.    for i := 1 to n do
  295.      begin
  296.       lo[i] := 1;
  297.       hi[i] := 1
  298.      end;
  299.  
  300.    order;
  301.    first;
  302.    numiter := 0;
  303.  
  304.    repeat
  305.      done := true;
  306.      numiter := succ(numiter);
  307.  
  308.      for i := 1 to n do
  309.            center[i] := 0.0;
  310.      for i := 1 to n do
  311.        if i<>hi[n] then
  312.         for j := 1 to m do
  313.               center[j] := center[j]+simp[i,j];
  314.  
  315.      for i := 1 to n do
  316.        begin
  317.          center[i] := center[i]/m;
  318.          next[i] := (1.0+alfa)*center[i]-alfa*simp[hi[n],i];
  319.        end;
  320.      sumresiduals(next);
  321.  
  322.      if next[n] <= simp[lo[n],n] then
  323.        begin
  324.          newvertex;
  325.          for i := 1 to m do
  326.            next[i] := gamma*simp[hi[n],i]+(1.0-gamma)*center[i];
  327.          sumresiduals(next);
  328.          if next[n] <= simp[lo[n],n] then newvertex;
  329.        end
  330.      else
  331.        begin
  332.          if next[n] <= simp[hi[n],n] then
  333.            newvertex
  334.          else
  335.            begin
  336.             for i := 1 to m do
  337.                 next[i] := beta*simp[hi[n],i]+(1.0-beta)*center[i];
  338.             sumresiduals(next);
  339.             if next[n] <= simp[hi[n],n] then
  340.               newvertex
  341.             else
  342.               begin
  343.                 for i := 1 to n do
  344.                   begin
  345.                     for j := 1 to m do
  346.                         simp[i,j] := (simp[i,j]+simp[lo[n],j])*beta;
  347.                   sumresiduals(simp[i]);
  348.                   end
  349.               end
  350.           end
  351.        end;
  352.        order;
  353.        for j := 1 to n do
  354.          begin
  355.            IF ABS(simp[hi[j], j]) > ABS(simp[lo[j], j]) {L.P.}
  356.            THEN {L.P.}
  357.                abs_max := ABS(simp[hi[j], j]) {L.P.}
  358.            ELSE {L.P.}
  359.                abs_max := ABS(simp[lo[j], j]); {L.P.}
  360.            error[j] :=
  361.                     (simp[hi[j],j]-simp[lo[j],j])/abs_max; {L.P.}
  362.            if done then
  363.            if error[j] > maxerr[j] then done := false;
  364.          end;
  365.  
  366.   until (done or (numiter = maxiter));
  367.  
  368.    for i := 1 to n do
  369.      begin
  370.        mean[i] := 0.0;
  371.        for j := 1 to n do
  372.             mean[i] := mean[i]+simp[j,i];
  373.        mean[i] := mean[i]/n;
  374.      end;
  375.   outputdata;
  376.   close(dout);
  377. { close(pltout); }
  378.   writeln;
  379.   writeln('Done....  results in file "SIMPLEX.OUT"');
  380. end.
  381. $